function [gesol,Phi_Nxz,prob_Nxz,Phi_NxYz,grid_Y,prob_NxYz,stats] = ...
    ComputeModelMoments_Model20220926_fast(parms,gesol_in)

    tic0 = tic;
    
    fprintf('ComputeModelMoments_Model20220926_fast: started computing moments.\n');

    gesol = gesol_in;
    
    %% steady state distributions

    % distribution of (N,L,z)
    [Phi_Nxz,prob_Nxz] = MakeSparseTransMatrix_Nxz(parms,gesol);
    Phi_Nxz_1y = mpower2(Phi_Nxz,12);
%     Phi_NLz_5y = Phi_NLz_1y^5;
    
    % check that the grid is sufficiently large
    prob_N = squeeze(sum(prob_Nxz,[2 3]));
    cdf_N = cumsum(prob_N(:));
    stats.grid_N_is_OK = cdf_N(2)<=0.001 && cdf_N(end-1)>=0.999;
    
    prob_x = squeeze(sum(prob_Nxz,[1 3]));
    cdf_x = cumsum(prob_x(:));
    stats.grid_x_is_OK = cdf_x(2)<=0.001 && cdf_x(end-1)>=0.999;
    
    % inflation process
    inflation_process.mu_pi = 6.1592e-04;
    inflation_process.rho_pi = 0.81;
    inflation_process.sigma_pi = 0.00245;
    inflation_process.theta_pi = -0.338;
    inflation_process.loading_dlogC = -0.035;
    gesol.inflation_process = inflation_process;
    
    % F(N,x,z,z') = -0.035(log C(N',x',z') - E[log C(N',x',z')|N,x,z])
    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;
    [NN_mesh,xx_mesh] = meshgrid(parms.grid_N(:),parms.grid_x(:));
    x_next_Nxzz = gen_law_of_motion_x(parms);
    Nnext_Nxz = gen_law_of_motion_N(parms,gesol);
    
    Pzz_Nxzz = permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]);
    C_next = zeros(NN,Nx,Nz,Nz); % C'(N,x,z,z')
    for iznext = 1:parms.Nz
        C_next(:,:,:,iznext) = reshape(interp2(NN_mesh,xx_mesh,gesol.C(:,:,iznext)',...
            Nnext_Nxz(:),x_next_Nxzz(:,iznext)),[NN Nx Nz]);
    end
    inflation_process.F_Nxzz = inflation_process.loading_dlogC*(C_next - repmat(sum(Pzz_Nxzz.*C_next,4),[1 1 1 Nz]));
    
    % distribution of (N,x,Y,z)
    NY = 5;
    [Phi_NxYz,grid_Y,prob_NxYz] = MakeSparseTransMatrix_NxYz(parms,gesol,inflation_process.rho_pi,...
        inflation_process.F_Nxzz,NY,prob_Nxz);
    
    grid_Y = grid_Y(:);
    
    prob_Y = squeeze(sum(prob_NxYz,[1 2 4]));
    cdf_Y = cumsum(prob_Y);
    stats.grid_Y_is_OK = cdf_Y(2)<=0.001 && cdf_Y(end-1)>=0.999;
    
    fprintf('Finished computing distributions. %g seconds elapsed.\n', toc(tic0));
    
    %% bond prices + entropy
    
    FLAG_CALC_ENTROPY = false;
    
    gesol = PriceRealBonds_Model20220926(parms,gesol,FLAG_CALC_ENTROPY);
    
    fprintf('Finished computing real bond prices. %g seconds elapsed.\n', toc(tic0));
    
    gesol.nom_bonds = PriceNomBonds_Model20220926(parms, gesol, inflation_process, 60);
    
    fprintf('Finished computing nominal bond prices. %g seconds elapsed.\n', toc(tic0));
    
%     % average entropy term structures
%     stats.entropy.horizon = gesol.entropy.horizon; 
%     NT = numel(gesol.entropy.horizon);
%     stats.entropy.MeanCondEntropy = ...
%         permute(reshape(gesol.entropy.CondEntropy,[NN*Nx*Nz, NT]),[2 1])*prob_Nxz(:);
%     MeanCondEntropy_given_z = zeros(NT,Nz);
%     for iz = 1:Nz
%         temp_entropy = permute(reshape(squeeze(gesol.entropy.CondEntropy(:,:,iz,:)),[NN*Nx, NT]),[2 1]);
%         temp_prob = prob_Nxz(:,:,iz);
%         MeanCondEntropy_given_z(:,iz) = temp_entropy*temp_prob(:)/sum(temp_prob(:));
%     end
%     stats.entropy.MeanCondEntropy_given_z = MeanCondEntropy_given_z;
%     
%     MeanCondEntropy_given_N = zeros(NT,NN);
%     for iN = 1:NN
%         temp_entropy = permute(reshape(squeeze(gesol.entropy.CondEntropy(iN,:,:,:)),[Nx*Nz, NT]),[2 1]);
%         temp_prob = squeeze(prob_Nxz(iN,:,:));
%         MeanCondEntropy_given_N(:,iN) = temp_entropy*temp_prob(:)/sum(temp_prob(:));
%     end
%     stats.entropy.MeanCondEntropy_given_N = MeanCondEntropy_given_N; % note this is only valid for regions of N that's visited!!
%     
%     MeanCondEntropy_given_x = zeros(NT,Nx);
%     for ix = 1:Nx
%         temp_entropy = permute(reshape(squeeze(gesol.entropy.CondEntropy(:,ix,:,:)),[NN*Nz, NT]),[2 1]);
%         temp_prob = squeeze(prob_Nxz(:,ix,:));
%         MeanCondEntropy_given_x(:,ix) = temp_entropy*temp_prob(:)/sum(temp_prob(:));
%     end
%     stats.entropy.MeanCondEntropy_given_x = MeanCondEntropy_given_x; % note this is only valid for regions of N that's visited!!
%     
%     % term structure of E[logM(t+H)-logM(t)]
%     MeanEdlogM_given_z = zeros(NT,Nz);
%     for iz = 1:Nz
%         temp_EdlogM = permute(reshape(squeeze(gesol.entropy.E_dlogM(:,:,iz,:)),[NN*Nx, NT]),[2 1]);
%         temp_prob = prob_Nxz(:,:,iz);
%         MeanEdlogM_given_z(:,iz) = temp_EdlogM*temp_prob(:)/sum(temp_prob(:));
%     end
%     stats.entropy.EdlogM_given_z = MeanEdlogM_given_z;
    
    %% stock price & returns
    gesol = PriceConsClaimUnlevered(parms,gesol,1e-7); % cum-dividend price
    
    fprintf('Finished computing stock prices. %g seconds elapsed.\n', toc(tic0));
    
%     NN = parms.NN; NL = parms.Nlambda; Nz = parms.Nz;
%     
    % stock returns
    FLAG_RUN_PARALLEL = true;
    stats = StockRetStats(parms,gesol,prob_Nxz,Phi_Nxz,stats,FLAG_RUN_PARALLEL);
    
    fprintf('Finished computing stock returns. %g seconds elapsed.\n', toc(tic0));
    
    %% macro moments
    
    U_Nxz = repmat(1 - parms.grid_N(:),[1 Nx Nz]);
    Y_Nxz = permute(repmat(exp(parms.grid_z(:)),[1 NN Nx]),[2 3 1]).*(1 - U_Nxz);

    [stats.U_mean, stats.U_std] = Calc_quarterly_ave_NHz(Phi_Nxz,prob_Nxz,U_Nxz);
    [stats.f_mean, stats.f_std] = Calc_quarterly_ave_NHz(Phi_Nxz,prob_Nxz,gesol.f);
    [stats.tightness_mean, stats.tightness_std,stats.tightness_AR1] = Calc_quarterly_AR1(Phi_Nxz,prob_Nxz,gesol.theta);

    stats.corr_U_V = Calc_quarterly_corr(Phi_Nxz,prob_Nxz(:),U_Nxz(:),gesol.theta(:).*U_Nxz(:));
    stats.corr_U_V_monthly = Calc_monthly_corr(prob_Nxz(:),U_Nxz(:),gesol.theta(:).*U_Nxz(:));

    stats.corr_U_Theta = Calc_quarterly_corr(Phi_Nxz,prob_Nxz(:),U_Nxz(:),gesol.theta(:));
    stats.corr_U_Theta_monthly = Calc_monthly_corr(prob_Nxz(:),U_Nxz(:),gesol.theta(:));

    [Y_mean, Y_std] = Calc_quarterly_ave_NHz(Phi_Nxz,prob_Nxz,Y_Nxz);
    [Nw_mean, Nw_std] = Calc_quarterly_ave_NHz(Phi_Nxz,prob_Nxz,gesol.wages.*(1-U_Nxz));
    stats.stdwagebill2stdY = Nw_std/Y_std;

    [stats.Cgrowth_AR1,stats.Cgrowth_vol] = moment_quarterly_growth_ar1(parms,gesol,gesol.C,Phi_Nxz,prob_Nxz);
    [stats.Ygrowth_AR1,stats.Ygrowth_vol] = moment_quarterly_growth_ar1(parms,gesol,Y_Nxz,Phi_Nxz,prob_Nxz);
    
    fprintf('Finished computing macro moments. %g seconds elapsed.\n', toc(tic0));
    
    
    %% cross sectional moments of consumption
    stats = CalcCrossSecMoments(parms, prob_Nxz, Phi_Nxz, stats);
    
    fprintf('Finished computing cross-sec moments. %g seconds elapsed.\n', toc(tic0));
    
    %% real yields
    mats = 12:12:60;
    real_yld_mean = zeros(size(mats));
    real_yld_std = zeros(size(mats));
%     real_yld_entropy_mean = zeros(size(mats));
%     real_yld_entropy_std = zeros(size(mats));
%     real_yld_EdlogM_mean = zeros(size(mats));
%     real_yld_EdlogM_std = zeros(size(mats));
    for n = 1:numel(mats)
        idx = find(gesol.real_bonds.maturities==mats(n));
        mat_yr = mats(n)/12;
        [real_yld_mean(n),real_yld_std(n)] = Calc_monthly_vol(prob_Nxz,-log(gesol.real_bonds.prices(:,:,:,idx))/mat_yr);
%         [real_yld_entropy_mean(n),real_yld_entropy_std(n)] = Calc_monthly_vol(prob_Nxz,-gesol.entropy.CondEntropy(:,:,:,idx)/mat_yr);
%         [real_yld_EdlogM_mean(n),real_yld_EdlogM_std(n)] = Calc_monthly_vol(prob_Nxz,-gesol.entropy.E_dlogM(:,:,:,idx)/mat_yr);
    end
    stats.real_yld_mean = 100*real_yld_mean;
    stats.real_yld_std = 100*real_yld_std;
%     stats.real_yld_entropy_mean = 100*real_yld_entropy_mean;
%     stats.real_yld_entropy_std = 100*real_yld_entropy_std;
%     stats.real_yld_EdlogM_mean = 100*real_yld_EdlogM_mean;
%     stats.real_yld_EdlogM_std = 100*real_yld_EdlogM_std;

    fprintf('Finished computing real yield moments. %g seconds elapsed.\n', toc(tic0));
    
    %% nominal yields
    nom_yld_mean = zeros(size(mats));
    nom_yld_std = zeros(size(mats));
    
    mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
    var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi + inflation_process.theta_pi^2)...
        *inflation_process.sigma_pi^2/(1 - inflation_process.rho_pi^2);
    var_e = inflation_process.sigma_pi^2;
    cov_Xe = inflation_process.sigma_pi^2;
    
% OLD VERSION WITHOUT PARFOR    
%     for n = 1:numel(mats)
%         idx = find(gesol.nom_bonds.maturities==mats(n));
%         mat_yr = mats(n)/12;
%         temp_input = permute(repmat(-log(gesol.nom_bonds.P_A(:,:,:,idx))/mat_yr,[1 1 1 NY]),[1 2 4 3]) ...
%             + permute(repmat(gesol.nom_bonds.P_C_a(idx)*grid_Y(:)/mat_yr,[1 NN Nx Nz]),[2 3 1 4]);
%         [temp_mean,temp_std] = Calc_monthly_vol(prob_NxYz,temp_input);
%         nom_yld_mean(n) = temp_mean + (gesol.nom_bonds.P_B_a(idx) + gesol.nom_bonds.P_B_b(idx)*mean_X)/mat_yr;
%         var_Xe_component = (gesol.nom_bonds.P_B_b(idx)^2*var_X + gesol.nom_bonds.P_B_c(idx)^2*var_e ...
%             + 2*gesol.nom_bonds.P_B_b(idx)*gesol.nom_bonds.P_B_c(idx)*cov_Xe)/(mat_yr^2);
%         nom_yld_std(n) = sqrt(temp_std^2 + var_Xe_component);
%     end

    % inputs into parfor
    P_A_list = zeros(NN, Nx, Nz, numel(mats));
    P_C_a_list = zeros(size(mats));
    P_B_a_list = zeros(size(mats));
    P_B_b_list = zeros(size(mats));
    P_B_c_list = zeros(size(mats));
    
    for n = 1:numel(mats)
        idx = find(gesol.nom_bonds.maturities==mats(n));
        P_A_list(:,:,:,n) = gesol.nom_bonds.P_A(:,:,:,idx);
        P_C_a_list(n) = gesol.nom_bonds.P_C_a(idx);
        P_B_a_list(n) = gesol.nom_bonds.P_B_a(idx);
        P_B_b_list(n) = gesol.nom_bonds.P_B_b(idx);
        P_B_c_list(n) = gesol.nom_bonds.P_B_c(idx);
    end

    parfor n = 1:numel(mats)
        mat_yr = mats(n)/12;
        temp_input = permute(repmat(-log(P_A_list(:,:,:,n))/mat_yr,[1 1 1 NY]),[1 2 4 3]) ...
            + permute(repmat(P_C_a_list(n)*grid_Y/mat_yr,[1 NN Nx Nz]),[2 3 1 4]);
        [temp_mean,temp_std] = Calc_monthly_vol(prob_NxYz,temp_input);
        nom_yld_mean(n) = temp_mean + (P_B_a_list(n) + P_B_b_list(n)*mean_X)/mat_yr;
        var_Xe_component = (P_B_b_list(n)^2*var_X + P_B_c_list(n)^2*var_e ...
            + 2*P_B_b_list(n)*P_B_c_list(n)*cov_Xe)/(mat_yr^2);
        nom_yld_std(n) = sqrt(temp_std^2 + var_Xe_component);
    end
    
    stats.nom_yld_mean = 100*nom_yld_mean;
    stats.nom_yld_std = 100*nom_yld_std;
    
    fprintf('Finished computing nominal yield moments. %g seconds elapsed.\n', toc(tic0));
    
    %% bond returns
    stats = BondRetStats(parms,gesol,prob_Nxz,Phi_Nxz_1y,stats);
    
    fprintf('Finished computing real bond stats. %g seconds elapsed.\n', toc(tic0));
    
    Phi_NxYz_1y = mpower2(Phi_NxYz,12);
    FLAG_RUN_PARALLEL = true;
    stats = BondRetStats_nominal(parms,gesol,prob_NxYz,Phi_NxYz_1y,grid_Y,stats, FLAG_RUN_PARALLEL);
    
    fprintf('Finished computing nominal bond stats. %g seconds elapsed.\n', toc(tic0));

    %% bond return decomposition into consumption, zeta, and cross terms

    % price fictitious bonds with M'=beta*(C'/C)^(-gamma)
    gesol_Mcons = gesol;
    gesol_Mcons.spd1period = gesol.spd1period_cons_comp;
    gesol_Mcons = PriceRealBonds_Model20220926(parms,gesol_Mcons);
    stats_Mcons = BondRetStats(parms,gesol_Mcons,prob_Nxz,Phi_Nxz_1y,[]);

    % price fictitious bonds with M'=zeta'/zeta
    gesol_Mzeta = gesol;
    gesol_Mzeta.spd1period = gesol.spd1period_zeta_comp;
    gesol_Mzeta = PriceRealBonds_Model20220926(parms,gesol_Mzeta);
    stats_Mzeta = BondRetStats(parms,gesol_Mzeta,prob_Nxz,Phi_Nxz_1y,[]);

    stats.bond_risk_prem_real_decomp_uncond.cons_term = stats_Mcons.bond_risk_prem_real_uncond;
    stats.bond_risk_prem_real_decomp_uncond.zeta_term = stats_Mzeta.bond_risk_prem_real_uncond;
    stats.bond_risk_prem_real_decomp_uncond.cross_term = stats.bond_risk_prem_real_uncond ...
        - stats_Mcons.bond_risk_prem_real_uncond - stats_Mzeta.bond_risk_prem_real_uncond;

    for iz = 1:parms.Nz
        stats.(['bond_risk_prem_real_decomp_z',num2str(iz)]).cons_term = ...
            stats_Mcons.(['bond_risk_prem_real_z',num2str(iz)]);
        stats.(['bond_risk_prem_real_decomp_z',num2str(iz)]).zeta_term = ...
            stats_Mzeta.(['bond_risk_prem_real_z',num2str(iz)]);
        stats.(['bond_risk_prem_real_decomp_z',num2str(iz)]).cross_term = ...
            stats.(['bond_risk_prem_real_z',num2str(iz)]) ...
            - stats_Mcons.(['bond_risk_prem_real_z',num2str(iz)]) ...
            - stats_Mzeta.(['bond_risk_prem_real_z',num2str(iz)]);
    end
    
    fprintf('Finished bond decomposition. %g seconds elapsed.\n', toc(tic0));
    
    %% predictive regressions, real & nominal
    
    FLAG_RUN_PARALLEL = true;
    [stats.betas_hpxr_theta_real, stats.betas_hpxr_fprob_real] = reg_hpxr_real(gesol,Phi_Nxz,prob_Nxz,FLAG_RUN_PARALLEL);
    
    fprintf('Finished predictive reg, real. %g seconds elapsed.\n', toc(tic0));
    
    % FLAG_RUN_PARALLEL = true; % process got killed
    FLAG_RUN_PARALLEL = false;
    [stats.betas_hpxr_theta_nom, stats.betas_hpxr_fprob_nom] = reg_hpxr_real_nom(gesol,Phi_NxYz,prob_NxYz,grid_Y,FLAG_RUN_PARALLEL);
    
    fprintf('Finished predictive reg, nominal. %g seconds elapsed.\n', toc(tic0));
    
    FLAG_RUN_PARALLEL = true;
    [stats.betas_slope_theta_real,stats.betas_yld_theta_real] = yield_regs_real(gesol, gesol.theta, prob_Nxz, FLAG_RUN_PARALLEL);
    
    fprintf('Finished yield reg, real. %g seconds elapsed.\n', toc(tic0));
    
    % FLAG_RUN_PARALLEL = true;
    FLAG_RUN_PARALLEL = false;
    [stats.betas_slope_theta_nom,stats.betas_yld_theta_nom] = yield_regs_nom(gesol, gesol.theta, grid_Y, prob_NxYz,FLAG_RUN_PARALLEL);
    
    fprintf('Finished yield reg, nominal. %g seconds elapsed.\n', toc(tic0));
    
    %% Fama-Bliss, real & nominal
    
    FLAG_RUN_PARALLEL = true;
    [stats.alphas_FB_real,stats.betas_FB_real,stats.R2s_FB_real, mats_yrs] = FBregs_real(gesol, Phi_Nxz, prob_Nxz, FLAG_RUN_PARALLEL);
    
    fprintf('Finished Fama-Bliss reg, real. %g seconds elapsed.\n', toc(tic0));
    
%     [stats.alphas_FB_real_ver2,stats.betas_FB_real_ver2,stats.R2s_FB_real_ver2, mats_yrs] = FBregs_real_ver2(gesol, Phi_Nxz, prob_Nxz);
    
    % FLAG_RUN_PARALLEL = true;
    FLAG_RUN_PARALLEL = false;
    [stats.alphas_FB_nom,stats.betas_FB_nom,stats.R2s_FB_nom, mats_yrs] = FBregs_nom(gesol, Phi_NxYz, prob_NxYz, grid_Y, FLAG_RUN_PARALLEL);
    
    fprintf('Finished Fama-Bliss reg, nominal. %g seconds elapsed.\n', toc(tic0));
    
%     [stats.alphas_FB_nom_ver2,stats.betas_FB_nom_ver2,stats.R2s_FB_nom_ver2, mats_yrs] = FBregs_nom_ver2(gesol, Phi_NxYz, prob_NxYz, grid_Y);

    %% Campbell-Shiller, real & nominal
%     [stats.alphas_CS_real,stats.betas_CS_real,stats.R2s_CS_real, mats_yrs] = CSregs_real(gesol, Phi_Nxz, prob_Nxz);
%     [stats.alphas_CS_nom,stats.betas_CS_nom,stats.R2s_CS_nom, mats_yrs] = CSregs_nom(gesol, Phi_NxYz, prob_NxYz, grid_Y);

    FLAG_RUN_PARALLEL = true;
    [stats.alphas_CS_real_ver2,stats.betas_CS_real_ver2,stats.R2s_CS_real_ver2, mats_yrs] = CSregs_ver2_real(gesol, Phi_Nxz, prob_Nxz, FLAG_RUN_PARALLEL);
    
    fprintf('Finished CS reg ver 2, real. %g seconds elapsed.\n', toc(tic0));
    
    % FLAG_RUN_PARALLEL = true;
    FLAG_RUN_PARALLEL = false;
    [stats.alphas_CS_nom_ver2,stats.betas_CS_nom_ver2,stats.R2s_CS_nom_ver2, mats_yrs] = CSregs_ver2_nom(gesol, Phi_NxYz, prob_NxYz, grid_Y, FLAG_RUN_PARALLEL);
    
    fprintf('Finished CS reg ver 2, nominal. %g seconds elapsed.\n', toc(tic0));

end